
%给定参数
Nt = 100;
p1 = [100, 100, 100, 100, 100, 100, 100, 100];
p2 = [50, 50, 50, 50, 50, 50, 50, 50];
p3 = [inf, inf, inf, inf, inf, inf, inf, inf];
h1 = 1000;
h2 = [h1/5, h1/3, h1/2, h1, h1*2, h1*3, h1*5, h1*24];
u2 = p2 ./ p1; u3 = p3 ./ p2;
v2 = h2 ./ h1;
i = sqrt(-1);
%进行计算
for ix = 1:8
    for  iy = 1: Nt
        B(ix,iy) = 10 ^ ((iy+10) / 15) / h1;
        C(ix,iy) = ( coth( 2*pi*sqrt(2) * exp( -i*pi/4) / B(ix,iy) + ...
            acoth( sqrt(u2(ix)) * coth( 2*pi*sqrt(2) * exp( -i*pi/4) / ...
            B(ix,iy) * v2(ix) / sqrt(u2(ix)) + acoth( sqrt( u3(ix)/u2(ix)))))))^2;
        D(ix,iy) = ( tanh( 2*pi*sqrt(2) * exp( -i*pi/4) / B(ix,iy) + ...
            atanh( sqrt(u2(ix)) * tanh( 2*pi*sqrt(2) * exp( -i*pi/4) / ...
            B(ix,iy) * v2(ix) / sqrt(u2(ix)) + atanh( sqrt( u3(ix)/u2(ix)))))))^2;
        Re(ix, iy) = real(C(ix, iy));          %实部
        Im(ix, iy) = imag(D(ix, iy));          %虚部
        A(ix, iy) = abs(C(ix, iy));
        E(ix, iy) = angle(C(ix,iy));
    end
    
end
%进行画图
% figure(9)
% loglog(B(1,:),Re1(1,:))
% hold on
% loglog(B(2,:),Re1(2,:))
% hold on
% loglog(B(3,:),Re1(3,:))
% hold on
% loglog(B(4,:),Re1(4,:))
% hold on
% loglog(B(5,:),Re1(5,:))
% hold on
% loglog(B(6,:),Re1(6,:))
% hold on
% loglog(B(7,:),Re1(7,:))
% hold on
% loglog(B(8,:),Re1(8,:))
% hold off
% title("三层大地电磁测深响应曲线 p3 = 100000")
% xlabel("lanbta1/h1"); ylabel("pw/p1")
% legend("1/5","1/3","1/2","1","2","3","5","24")

figure(9)
loglog(B(1,:),A(1,:))
hold on
loglog(B(2,:),A(2,:))
hold on
loglog(B(3,:),A(3,:))
hold on
loglog(B(4,:),A(4,:))
hold on
loglog(B(5,:),A(5,:))
hold on
loglog(B(6,:),A(6,:))
hold on
loglog(B(7,:),A(7,:))
hold on
loglog(B(8,:),A(8,:))
hold off
title("三层大地电磁测深响应曲线 p3 = 100000")
xlabel("lanbta1/h1"); ylabel("pw/p1")
legend("1/5","1/3","1/2","1","2","3","5","24")
figure(10)
loglog(B(1,:),E(1,:))
hold on
loglog(B(2,:),E(2,:))
hold on
loglog(B(3,:),E(3,:))
hold on
loglog(B(4,:),E(4,:))
hold on
loglog(B(5,:),E(5,:))
hold on
loglog(B(6,:),E(6,:))
hold on
loglog(B(7,:),E(7,:))
hold on
loglog(B(8,:),E(8,:))
hold off
